#20210607-CK与T一一对应
#20210525-分位数计算CI


rm(list=ls())

library(xlsx)
library(plyr)
library(RColorBrewer)
library(ggplot2)
library(scales)
library(dplyr)
library(ggpubr)
library(gtable)
library(grid)
library(ggExtra)
library(cowplot)

######################################################
###     数据读取与设置  Data loading and setting   ###
######################################################
Data <- read.xlsx("E:\\R\\working_area\\data.xlsx",sheetName = "data1",encoding = "UTF-8")  #载入数据Load data
Data$PAIR[Data$PAIR == "A_N"] <- "Native"
Data$PAIR[Data$PAIR == "N_A"] <- "Invasive"

n_iter <- 1000 #迭代次数 iteration
n_samp <- 100 #抽样次数 Number of samples

t_gcf <- as.vector(unique(Data$F)) #全球变化因子类型
t_index <- c("Performance","Competitiveness") #指标类型
t_plant <- c("Native","Invasive") #植物类型

type <- plyr::rename(expand.grid(t_gcf, t_index, t_plant), c(Var1 = "gcf", Var2 = "index", Var3 = "plant")) #排列组合
#type$remake <- paste0(type$gcf,'_', type$index, '_', type$plant) #汇总
type$n_gcf <- nchar(as.vector(type$gcf))


#######################################
###     数据处理  data processing   ###
#######################################

#文献与研究  Literature and research
literature <- list()
literature[[1]] <- read.xlsx("E:\\R\\working_area\\data.xlsx",sheetName = "literature",colIndex = c(1:5)) #文献国别 Literature country & 研究国别 Research country
literature[[2]] <- read.xlsx("E:\\R\\working_area\\data.xlsx",sheetName = "literature",colIndex = c(6:8))  #文献、研究年份 Literature & Research year
literature[[3]] <- read.xlsx("E:\\R\\working_area\\data.xlsx",sheetName = "literature",colIndex = c(12:16)) #年份、国别 Literature & Research year & country
literature[[4]] <- read.xlsx("E:\\R\\working_area\\data.xlsx",sheetName = "literature",colIndex = c(12:14))  #文献、研究类型 Literature & Research type


#计算样本量 Sample size
N <- rbind(plyr::rename(ddply(Data,c("F","PAIR"),summarize,N = length(F)),c(F="gcf",PAIR="plant")),plyr::rename(ddply(Data,c("Type_F","PAIR"),summarize,N = length(F)),c(Type_F="gcf",PAIR="plant")))
N$gcf[N$gcf=="1"] <- "GCF(S)"
N$gcf[N$gcf=="2"] <- "GCF(M)"


#拔靴法计算真实数据 Actual data from bootStrap
actual_data <- list()
for (i in c(1:8)) {
  actual_data[[nrow(type)+i]] <- data.frame()
}
for (i in c(1:nrow(type))) {
  sample <- data.frame(plant = Data$PAIR[Data$PAIR == type$plant[i] & Data$F == type$gcf[i]],
                       gcf = Data$F[Data$PAIR == type$plant[i] & Data$F == type$gcf[i]],
                       index = type$index[i])
  
  
  if(type$index[i] == "Performance") {
    sample$TR <- Data$GT[Data$PAIR == type$plant[i] & Data$F == type$gcf[i]]
    CT <- Data$NT[Data$PAIR == type$plant[i] & Data$F == type$gcf[i]]
  }
  if(type$index[i] == "Competitiveness") {
    sample$TR <- Data$GR[Data$PAIR == type$plant[i] & Data$F == type$gcf[i]]
    CT <- Data$NR[Data$PAIR == type$plant[i] & Data$F == type$gcf[i]]
  }
  bs <- c()
  for (j in c(1:n_iter)) {
  set.seed(j)
  bs <- append(bs, mean(sample(sample$TR, n_samp, replace = TRUE), na.rm = T) - mean(sample(CT, n_samp, replace = TRUE), na.rm = T))
}
  actual_data[[i]] <- data.frame(plant = rep(sample[1,1], n_iter), gcf = rep(sample[1,2], n_iter), index = rep(sample[1,3], n_iter), effect_size = bs)
  if (type$index[i] == "Performance" & type$plant[i] == "Native" & type$n_gcf[i] == 1) {
    actual_data[[nrow(type)+1]] <- rbind(actual_data[[nrow(type)+1]], actual_data[[i]])
    actual_data[[nrow(type)+1]]$gcf <- "GCF(S)"
  }
  if (type$index[i] == "Performance" & type$plant[i] == "Native" & type$n_gcf[i] == 2) {
    actual_data[[nrow(type)+2]] <- rbind(actual_data[[nrow(type)+2]], actual_data[[i]])
    actual_data[[nrow(type)+2]]$gcf <- "GCF(M)"
  }
  if (type$index[i] == "Performance" & type$plant[i] == "Invasive" & type$n_gcf[i] == 1) {
    actual_data[[nrow(type)+3]] <- rbind(actual_data[[nrow(type)+3]], actual_data[[i]])
    actual_data[[nrow(type)+3]]$gcf <- "GCF(S)"
  }
  if (type$index[i] == "Performance" & type$plant[i] == "Invasive" & type$n_gcf[i] == 2) {
    actual_data[[nrow(type)+4]] <- rbind(actual_data[[nrow(type)+4]], actual_data[[i]])
    actual_data[[nrow(type)+4]]$gcf <- "GCF(M)"
  }  
  if (type$index[i] == "Competitiveness" & type$plant[i] == "Native" & type$n_gcf[i] == 1) {
    actual_data[[nrow(type)+5]] <- rbind(actual_data[[nrow(type)+5]], actual_data[[i]])
    actual_data[[nrow(type)+5]]$gcf <- "GCF(S)"
  }
  if (type$index[i] == "Competitiveness" & type$plant[i] == "Native" & type$n_gcf[i] == 2) {
    actual_data[[nrow(type)+6]] <- rbind(actual_data[[nrow(type)+6]], actual_data[[i]])
    actual_data[[nrow(type)+6]]$gcf <- "GCF(M)"
  }  
  if (type$index[i] == "Competitiveness" & type$plant[i] == "Invasive" & type$n_gcf[i] == 1) {
    actual_data[[nrow(type)+7]] <- rbind(actual_data[[nrow(type)+7]], actual_data[[i]])
    actual_data[[nrow(type)+7]]$gcf <- "GCF(S)"
  }
  if (type$index[i] == "Competitiveness" & type$plant[i] == "Invasive" & type$n_gcf[i] == 2) {
    actual_data[[nrow(type)+8]] <- rbind(actual_data[[nrow(type)+8]], actual_data[[i]])
    actual_data[[nrow(type)+8]]$gcf <- "GCF(M)"
  }
}
rm(i, j, bs, sample, CT)


#由真实值估算 Calculated data from actual data
calculated_data <- list()
for (i in c(1:4)) {
  calculated_data[[nrow(type)+i]] <- data.frame()
}
for (i in which(type$n_gcf == 2)) {
  gcf1 <- which(type$gcf == substr(as.vector(type$gcf[i]), 1, 1) &
                  type$index == as.vector(type$index[i]) &
                  type$plant == as.vector(type$plant[i]))
  gcf2 <- which(type$gcf == substr(as.vector(type$gcf[i]), 2, 2) &
                  type$index == as.vector(type$index[i]) &
                  type$plant == as.vector(type$plant[i]))
  calculated_data[[i]] <- cbind(actual_data[[i]][-4],
                                effect_size = actual_data[[gcf1]]$effect_size + actual_data[[gcf2]]$effect_size)
  calculated_data[[i]] <- cbind(calculated_data[[i]],
                                difference_size = actual_data[[i]]$effect_size - calculated_data[[i]]$effect_size)
  if (type$index[i] == "Performance" & type$plant[i] == "Native") {
    calculated_data[[nrow(type)+1]] <- rbind(calculated_data[[nrow(type)+1]], calculated_data[[i]])
    calculated_data[[nrow(type)+1]]$gcf <- "GCF(M)"
  }
  if (type$index[i] == "Competitiveness" & type$plant[i] == "Native") {
    calculated_data[[nrow(type)+2]] <- rbind(calculated_data[[nrow(type)+2]], calculated_data[[i]])
    calculated_data[[nrow(type)+2]]$gcf <- "GCF(M)"
  }
  if (type$index[i] == "Performance" & type$plant[i] == "Invasive") {
    calculated_data[[nrow(type)+3]] <- rbind(calculated_data[[nrow(type)+3]], calculated_data[[i]])
    calculated_data[[nrow(type)+3]]$gcf <- "GCF(M)"
  }
  if (type$index[i] == "Competitiveness" & type$plant[i] == "Invasive") {
    calculated_data[[nrow(type)+4]] <- rbind(calculated_data[[nrow(type)+4]], calculated_data[[i]])
    calculated_data[[nrow(type)+4]]$gcf <- "GCF(M)"
  }
}
calculated_data <- calculated_data[-which(sapply(calculated_data, is.null))]
rm(i, gcf1, gcf2)

#计算合并值和CI以及过“0”概率
result_eff <- data.frame(plant = NA, gcf = NA, index = NA, mean = NA, down = NA, up = NA, p_value = NA)
for (i in c(1:length(actual_data))) {
  result_eff[i,c(1:3)] <- actual_data[[i]][1,-4]
  result_eff$mean[i] <- mean(actual_data[[i]]$effect_size)
  result_eff$down[i] <- quantile(actual_data[[i]]$effect_size, .025)
  result_eff$up[i] <- quantile(actual_data[[i]]$effect_size, .975)
  result_eff$p_value[i] <- length(which(actual_data[[i]]$effect_size>0))/length(actual_data[[i]]$effect_size)
}
result_eff <- merge(result_eff,N,by = c("gcf","plant"), all=TRUE)
result_eff$index[result_eff$index == 1] <- "Performance"
result_eff$index[result_eff$index == 2] <- "Competitiveness"
result_eff$index <- factor(result_eff$index, levels=c("Performance", "Competitiveness"))
result_eff$N[result_eff$down/result_eff$up > 0] <- paste0(result_eff$N[result_eff$down/result_eff$up > 0], " *")


result_diff <- data.frame(plant = NA, gcf = NA, index = NA, mean = NA, down = NA, up = NA, p_value = NA)
for (i in c(1:length(calculated_data))) {
  result_diff[i,c(1:3)] <- calculated_data[[i]][1,c(-4,-5)]
  result_diff$mean[i] <- mean(calculated_data[[i]]$difference_size)
  result_diff$down[i] <- quantile(calculated_data[[i]]$difference_size, .025)
  result_diff$up[i] <- quantile(calculated_data[[i]]$difference_size, .975)
  result_diff$p_value[i] <- length(which(calculated_data[[i]]$difference_size>0))/length(calculated_data[[i]]$difference_size)
}
result_diff$index[result_diff$index == 1] <- "Performance"
result_diff$index[result_diff$index == 2] <- "Competitiveness"
result_diff$index <- factor(result_diff$index, levels=c("Performance", "Competitiveness"))
result_diff$N[result_diff$down/result_diff$up > 0] <- "*"

write.xlsx(result_eff,"E:\\R\\working_area\\result.xlsx",sheetName = "eff")
write.xlsx(result_diff,"E:\\R\\working_area\\result.xlsx",sheetName = "dif",append = TRUE)


##############################
###     绘图  plotting     ###
##############################

windowsFonts(myFont = windowsFont("Times New Roman")) #字体修改

# Compute the position of labels
F0 <- list()
F0[[11]] <- ggplot(literature[[1]]) +
  geom_col(aes(x=Code, y=Literature)) +
  coord_flip() +
  theme_void()

F0[[21]] <- ggplot(literature[[1]]) +
  geom_col(aes(x=Code, y=Study)) +
  coord_flip() +
  theme_void()

F0[[12]] <- ggplot(literature[[2]]) +
  geom_col(aes(x=Year, y=Literature)) +
  theme_void()

F0[[22]] <- ggplot(literature[[2]]) +
  geom_col(aes(x=Year, y=Study)) +
  theme_void()

F0[[13]] <- ggplot(literature[[3]][literature[[3]]$Group == "Literature",], aes(x=Year, y=Code, size=Number)) +
  geom_point() +
  theme_bw() +
  theme(legend.position="bottom") +
  scale_y_continuous(breaks=literature[[3]]$Code, labels = literature[[3]]$Country) +
  scale_x_continuous(limits=c(2005,2019),breaks = c(2005,2007,2009,2011,2013,2015,2017,2019)) +
  labs(x = "", y= "") +
  theme(title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12))

F0[[23]] <- ggplot(literature[[3]][literature[[3]]$Group == "Study",], aes(x=Year, y=Code, size=Number)) +
  geom_point() +
  theme_bw() +
  theme(legend.position="bottom") +
  scale_y_continuous(breaks=literature[[3]]$Code, labels = literature[[3]]$Country) +
  scale_x_continuous(limits=c(2005,2019),breaks = c(2005,2007,2009,2011,2013,2015,2017,2019))+
  labs(x = "", y= "") +
  theme(title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12))


F0[[1]] <- cowplot::ggdraw() + 
  cowplot::draw_plot(F0[[12]], 0.18, 0.8, 0.62, 0.2) + 
  cowplot::draw_plot(F0[[13]], 0, 0, 0.8, 0.8) +
  cowplot::draw_plot(F0[[11]], 0.8, 0.2, 0.2, 0.6) 

F0[[2]] <- cowplot::ggdraw() + 
  cowplot::draw_plot(F0[[22]], 0.18, 0.8, 0.62, 0.2) + 
  cowplot::draw_plot(F0[[23]], 0, 0, 0.8, 0.8) +
  cowplot::draw_plot(F0[[21]], 0.8, 0.2, 0.2, 0.6)



F1 <- ggplot(data=result_eff[result_eff$gcf %in% c("GCF(S)","C","N","D","P","T"),])+
  aes(x = mean, y = gcf, color = plant)+
  geom_errorbarh(aes(xmin = down, xmax = up), height = 0.2, size = 0.7, position = position_dodge(width = -0.5))+
  geom_point(size = 3, position = position_dodge(width = -0.5), aes(shape = plant))+
  geom_vline(xintercept = 0, linetype = "dashed", size = 0.7)+
  scale_color_brewer(palette="Set1")+
  geom_text(data=result_eff[result_eff$gcf %in% c("GCF(S)","C","N","D","P","T"),],aes(x = up+0.08, y = gcf, color = plant,label = N),position=position_dodge(width = -0.5))+
  theme_bw()+
  theme(axis.title.y=element_blank())+
  theme(legend.title=element_blank(),legend.text=element_text(family = "myFont", size = 14),legend.justification=c(0.01,0.01), legend.position=c(0.01,0.01))+
  theme(title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12))+
  scale_y_discrete(limits = factor(c("GCF(S)","C","N","D","P","T"), levels=c("GCF(S)","C","N","D","P","T"))) +
  labs(x = "Effect size")+
  #scale_x_continuous(limits=c(-1.1,1.1),breaks = c(-0.75,0,0.75))+
  facet_wrap(vars(index))+
  theme(strip.text.x = element_text(family = "myFont",size=14))

F2 <- ggplot(data=result_eff[result_eff$gcf %in% c("GCF(M)","NC","DC","DN","PN","TC","TN","TD"),])+
  aes(x = mean, y = gcf, color = plant)+
  geom_errorbarh(aes(xmin = down, xmax = up), height=0.5, size=0.7 ,position = position_dodge(width = -0.5))+
  geom_point(size = 3, position = position_dodge(width = -0.5), aes(shape = plant))+
  geom_vline(xintercept = 0, linetype = "dashed", size= 0.7)+
  scale_color_brewer(palette="Set1")+
  geom_text(data=result_eff[result_eff$gcf %in% c("GCF(M)","NC","DC","DN","PN","TC","TN","TD"),],aes(x=up+0.2,y=gcf, color = plant,label=N),position=position_dodge(width = -0.5))+
  theme_bw()+
  theme(axis.title.y=element_blank())+
  theme(legend.title=element_blank(),legend.text=element_text(family = "myFont", size = 14),legend.justification=c(0.01,0.01), legend.position=c(0.01,0.01))+
  theme(title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12))+
  scale_y_discrete(limits = factor(c("GCF(M)","NC","DC","DN","PN","TC","TN","TD"), levels=c("GCF(M)","NC","DC","DN","PN","TC","TN","TD"))) +
  labs(x = "Effect size")+
  #scale_x_continuous(limits=c(-1.7,0.8),breaks = c(-1.5,-1.0,-0.5,0,0.5))+
  facet_wrap(vars(index))+
  theme(strip.text.x = element_text(family = "myFont",size=14))

F3 <- ggplot(data=result_diff[result_diff$gcf %in% c("GCF(M)","NC","DC","DN","PN","TC","TN","TD") ,])+
  aes(x = mean, y = gcf, color = plant)+
  geom_errorbarh(aes(xmin = down, xmax = up), height=0.5, size=0.7 ,position = position_dodge(width = -0.5))+
  geom_point(size = 3, position = position_dodge(width = -0.5), aes(shape = plant))+
  geom_vline(xintercept = 0, linetype = "dashed", size= 0.7)+
  scale_color_brewer(palette="Set1")+
  geom_text(data=result_diff[result_diff$gcf %in% c("GCF(M)","NC","DC","DN","PN","TC","TN","TD") ,], aes(x = up+0.1, y = gcf, color = plant,label = N),position=position_dodge(width = -0.5))+
  theme_bw()+
  theme(axis.title.y=element_blank())+
  theme(legend.title=element_blank(),legend.text=element_text(family = "myFont", size = 14),legend.justification=c(0.01,0.01), legend.position=c(0.01,0.01))+
  theme(title=element_text(family = "myFont",size=14,),
        axis.text.x = element_text(family = "myFont",size=12),
        axis.text.y = element_text(family = "myFont",size=12))+
  scale_y_discrete(limits = factor(c("GCF(M)","NC","DC","DN","PN","TC","TN","TD"), levels=c("GCF(M)","NC","DC","DN","PN","TC","TN","TD"))) +
  labs(x = "Difference size")+
  scale_x_continuous(limits=c(-1.5,1.2), breaks = c(-1.0,-0.5,0,0.5,1.0))+
  facet_wrap(vars(index))+
  theme(strip.text.x = element_text(family = "myFont",size=14))

jpeg(file = "Fig_3.jpg",width = 2100,height = 2700,units = "px",res = 300,quality = 100)
print(ggpubr::ggarrange(F0[[1]],F0[[2]], labels = c("Literature", "Study"),nrow = 2))
dev.off()

jpeg(file = "Fig_4.jpg",width = 2400,height = 1600,units = "px",res = 300,quality = 100)
F1
dev.off()

jpeg(file = "Fig_5.jpg",width = 2400,height = 1600,units = "px",res = 300,quality = 100)
F2
dev.off()

jpeg(file = "Fig_6.jpg",width = 2400,height = 1600,units = "px",res = 300,quality = 100)
F3
dev.off()


